Collapse and revival of excitations in Bose-Einstein condensates 
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We study the energies and decay of elementary excitations in weakly interacting Bose-Einstein 
condensates within a finite temperature gapless second order theory. The energy shifts for the 
high-lying collective modes turn out to be systematically negative compared with the Hartree- 
Fock-Bogoliubov-Popov approximation and the decay of the low-lying modes are found to exhibit 
collapse and revival effects. In addition, perturbation theory is used to qualitatively explain the 
experimentally observed Beliaev decay process of the scissors mode. 
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I. INTRODUCTION 

The partially Bose-Einstein condensed trapped atomic 
gases provide an excellent testbench for developing finite 
temperature quantum theories. These weakly interacting 
systems can be modelled from first principles, and the 
experiments yield accurate and detailed information for 
comparison. Especially, the energies and decay rates of 
low-energy collective excitations have been measured at 
different temperatures and the results provide stringent 
tests for theoretical models. 

For dilute condensates at temperatures much lower 
than the condensation temperature T^, the Bogoliubov 
approximation consisting of the Gross-Pitaevskii (GP) 
equation for the condensate wave function and the Bo- 
goliubov equations for the quasiparticle excitations has 
proven to be accurate in describing the collective modes 
of the system. For higher temperatures one has to 
take into account the effects of the thermal gas com- 
ponent. Developing a theory that is computationally 
feasible and correctly models the system at tempera- 
tures approaching Tc is a challenging task. The most 
commonly used finite-temperature theory is the Hartree- 
Fock-Bogoliubov-Popov (HFB-Popov) approximation. It 
neglects the dynamics of the thermal gas and the modi- 
fications in particle correlations induced by the conden- 
sate, but predicts quasiparticle energies in fair agreement 
with the experiments [1]. The energy of the quadrupole 
modes having azimuthal angular momentum quantum 
numbers qg — ±2 deviates from the theoretical predic- 
tion for temperatures above 0.6Tc, but lately this devia- 
tion has been interpreted to mainly arise from improper 
modelling of the time-dependent external potential used 
in the experiments to excite the collective modes [2]. 

In order to take into account the leading order quasi- 
particle interactions and the correlations induced by the 
condensate in the inhomogeneous case, several theoreti- 
cal approaches have been suggested [3-10]. The dynam- 
ics of the condensate and the thermal gas has also been 
studied using various kinetic theories [11-15]. The sec- 
ond order theory for inhomogeneous, partially condensed 
gases presented in Refs. [9, 10] uses systematic perturba- 
tion theory to take into account the interaction terms 
in the Hamiltonian. Recently, this theory was extended 



to take into account the time-dependent external per- 
turbation used to drive the system in the experiments, 
leading to an agreement with the measured energies and 
the damping rates of the collective modes [2, 16]. 

The second order theories are computationally chal- 
lenging, and there has been only a few numerical investi- 
gations of their predictions [2, 17, 18]. In this paper, we 
calculate the spectral distributions of the quasiparticle 
energies for a partially condensed Bose-Einstein conden- 
sate (BEG), and compare the quasiparticle energies to 
the HFB-Popov results as functions of temperature. Es- 
pecially, we analyze the quasiparticle dynamics implied 
by the spectral distributions, observing that some col- 
lective modes should exhibit notable collapse and revival 
effects in trapped condensates. The possible existence 
of this phenomenon has been pointed out previously in 
Refs. [9, 10] (see also Ref. [19]), but, however, it has not 
been studied in detail before. The collapse and revival of 
the excitations indicates that the energies and the damp- 
ing rates alone do not suffice to describe the dynamics of 
these modes, i.e., the commonly used damped sinusoidal 
fit to the experimental data may not be sufficient to de- 
scribe the longer term dynamics of some modes. 

The structure of the paper is the following: In Section 
II, we describe the second order theory used in the analy- 
sis. Section III is devoted to a discussion of the numerical 
methods used to calculate the excitation spectra and the 
dynamics of modes. In Section IV, we analyze the second 
order corrections to the excitation spectrum as functions 
of temperature, and in Section V we study the decay 
of certain modes. Section VI consists of discussion and 
summary of the results. 



II. SECOND ORDER THEORY 

In this section, following Refs. [9, 10], we present the 
second order formalism for calculating the quasiparticle 
spectral distributions for a partially condensed, dilute, 
trapped BEG at finite temperatures. The starting point 
is the usual second quantized Hamiltonian for structure- 
less bosons 



H 



ijkm 



(ij\V\km)ala]dkam, (1) 
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where the creation and annihilation operators for a par- 
ticle in state \i) are denoted by a| and a^, respectively. 
The single particle Hamiltonian is given by the sum of 
the kinetic energy and the external trapping potential as 



2m 



Vi;rap(r), 



and the dominant s-wave scattering at low temperatures 
can be modelled by the effective low energy interaction 
potential 



V{v) = 



(2) 



where a is the scattering length and m the atomic mass. 
This effective potential is inapplicable at high energies 
and leads to ultraviolet divergences in the theory which 
have to be renormalized in a proper way, see Appendix. 

We choose to use a canonical ensemble with fixed total 
number of particles A'^. By defining the bosonic number 
conserving operators dj = [{Nq + l)^^/-^do]^aj, where the 
index refers to the condensate state and Nq = ajao, 
one can write the Hamiltonian (1) as 



H = ^Hi + OiNo[S/Nof'), 



i=0 



where 
Ho 
Hi = 

Ho = 



No 



(o|iJ"p|o) + -iVo(oo|y''|oo) 



(3) 



(4) 



^oJ2[{i\H'P\0)+No{iO\V'\00)\ dt+h.c, (5) 



^ [{i\H'P\j) - XSij + 2iVo(0i|y^|j0) 



No,. 



(uiv^ioo), 



i,tAt 



h.C. 



H3=J2 [VN'o{ij\V''\kO)ala]ak+h.c. 



H4= ^ ^{ij\V''\km)alalakar, 



A(iVex), (6) 
(7) 
(8) 



and S = Ncx — (A^cx) is the number fluctuation operator 
of the noncondensate particles. The symmetrized matrix 
elements of the two-particle interaction potential V{r) 
are defined as 

{ij\V'\km) = h{ij\V\km) + {ji\V\km)], 



and A as 



A = (OliJ^'PlO) +7Vo(00|F''|00), 



where the average number of atoms in the condensate 
state is given hy No = N — (iVex)- Above the averages 



(...) refer to quantum expectation values and h.c. stands 
for hermitian conjugate. 

In the zcroth order approximation, one solves the 
ground state |0) of Ho alone, which makes the linear 
Hamiltonian Hi to vanish. The excitations are found 
in lowest order by diagonalizing H2 and the number of 
the condensed particles A^o has to be tuned such that the 
total number of particles satisfies N = No + A^ex- 

It is convenient the use an orthonormal single-particle 
basis ^i(r) = (r|i) for alH = 0, 1, . . . , where Co(r) is the 
condensate wave function given by the Gross-Pitaevskii 
equation 

- 2^V2Co + Ft,ap(r)Co + A^oC/olCoWl'Co = ACo(r), (9) 

with Uq = 'inah^ /m. The GP equation is obtained 
by minimizing (Ho) with respect to Co(r). Diagonal- 
izing H2 using the Bogoliubov transformation pi = 



tions 



£(r) M{r) 
-M*{r) -C{r) 



results in the Bogoliubov equa- 



Ui{r) 



Ui{r) 



(10) 



where Ui{r) = Y,.^^ UijCj{r) and Vi{r) = J2jjio VijQ{r) 
are the quasiparticle amplitudes, e, the quasiparticle en- 
ergies, and operators C{r) = ff'^P — A + 2A"ot^o|Co(r)P and 
A^(r) = A^oC^oCo('') have been introduced. The quasi- 
particle amplitudes must satisfy the orthogonality and 
symmetry relations 



/ 
/ 



dr[ui{r)u*{r) - Vi{r)Vj{r)] = 6ij, 
dr[ui{r)vj{r) -Vi{r)uj{r)] = 0, 



(11) 



for the Bogoliubov transformation to be canonical. The 
quasiparticles must also be orthogonal to the conden- 
sate state, i.e., / drQ{r)ui{r) = / drCo(r)wi(r) = 0. 
The Bogoliubov equations have the zero-energy solution 
{uo{r),Vo{r)} = {Co{r), -Q{r)}, and projection to this 
homogeneous solution should always be subtracted from 
the quasiparticle amplitudes. 

To calculate the next lowest-order mean fields, 
i.e., the density of the thermal atoms p(r) = 
Sij^^o Cj (r)Ci(r)(dj(3!i) and the so-called anomalous av- 
erage K{r) = J2ijjto 0(r)Cj(r)(djdj), we express the par- 
ticle operators in terms of the quasiparticle operators, 
yielding 

p(r) = ^{[\u,{r)f + \vi{r)\']ni + \vi{r)f}, (12) 
K{r) = J2MrK{r){2ni + l). (13) 

In principle, the quasiparticle populations n, = 0l(3i) 
should be calculated from the requirement that the 
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canonical partition function — X^jn;} ^ /^^'({"'i) min- 
imizes the free energy T = —k^T log Zc- However, to a 
good approximation [20] one may use the non-interacting 
quasiparticle gas result rii = {z~^e^^^ — where the 

fugacity is calculated from the relation z = Nq/{1 + Nq). 

In calculating the perturbative corrections to the 
zeroth-order theory corresponding to Eqs. (9) and (10), 
it is convenient to first calculate the improved condensate 
wave frmction Co (r) from the generalized Gross-Pitaevskii 
(GGP) equation 

- I^V^CoW + 14rap(r)Co(r) + 7VoC/o|Co(r)|2Co(r) 
+ 2[/o/)(r)Co(r) + ?7o/«(r)Co*(r) = AgCo(r), (14) 

which is obtained by minimizing (Hq) + {H^)- Expressing 
the terms in the Hamiltonian as 

Hi = Hi[Co]+^Hi, (15) 

where 

/^H, = H,[Qo\ - HQo], (16) 
one finds the perturbative Hamiltonian 

^pert = AHo + Ai?i + AJf2 +Hz+ Hi, (17) 

where the non-quadratic terms and H4 are to be cal- 
culated using the improved condensate wave function. 
Note that our notation for AHi differs somewhat from 
that in Refs. [9, 10]. 

The perturbation term AHq is just a real number and 
can be easily taken into account. In addition to it, in 
first order perturbation theory only the terms AH2 and 
H4 containing even numbers of quasiparticle operators 
contribute to the energy shift 

Epertis, 1) = (sl^pertjs), (18) 

where \s) is a quasiparticle occupation number eigen- 
state. In second order perturbation theory, one can in 
fact neglect the terms AH2 and H4, because it turns out 
that their contribution is of the same order as the con- 
tribution of the other terms in third order perturbation 
theory [9, 10]. Thus, one only needs to calculate 

The quasiparticle energies are calculated as total en- 
ergy changes in the system when the corresponding quasi- 
particle occupation number is increased by one, while the 
total number of particles is held constant. This yields the 
corrected excitation energy 

E,{z') = e, + AEl + Ai?f,^p^ + AEl + Ai?f (/), (20) 

where the A-terms are given in Eqs. (A.3)-(A.6) and the 
complex energy parameter z' should not be mixed with 



the fugacity. Calculating the excitation energies as func- 
tions of z' yields the dynamics of the excitations in the 
following way: The time evolution operator U{t) of the 
system may be written in terms of the Fourier transform 
of the resolvent operator G{z') = {z' — H)~^ as [21] 

U{t) = / e-^'^*Im[G(?iw - iQ)\dLO. (21) 

7^ J -00 

Let us define^ the projection of the resolvent to state p as 
Gp{z) = {p\G\p), which may be approximated to second 
order as 

Gpioj) = [huj - Ep{Tvuj)\-'^. (22) 

Finally, it is seen that the imaginary part of the projected 
resolvent Fp{u} — zO) = Im[Gp(w — iO)] gives the spectral 
distribution of the mode p and the Fourier transform of 
-Pp(w) yields its time dependence. 

The need to calculate quasiparticle energies as func- 
tions of z' is naturally related to the fact that one takes 
into account quasiparticle interactions, though only to 
the lowest order, and the quasiparticle states are no more 
energy eigcnstates having infinite lifetime. In addition, 
in computing the quasiparticle energies z' must have 
a small imaginary part acting as a regularizer for the 
otherwise divergent expressions for second order energy 
shifts. One may note that setting z' = Cp yields the 
usual Rayleigh-Schrodinger perturbation theory, while 
the Brillouin-Wigner perturbation theory corresponds to 
solving the equation Ep{z') = z' . 

In conclusion, the second order theory may be used to 
calculate the energies and the dynamics of quasiparticles. 
First the GP equation (9) is solved together with the Bo- 
goliubov equations (10) for a given total particle number 
N = No + Nc-y,. Then the GGP equation (14) is solved, 
after which the spectrum -Fp(w) may be extracted for 
each excitation p using the energy corrections presented 
in Eqs. (A.3)-(A.6). In addition, one has to take care of 
proper ultraviolet renormalization; The quantities k{t) 
and AE^ are to be replaced by their renormalized values 
given in Eqs. (A. 2) and (A. 9) in all calculations. 

III. NUMERICAL METHODS 

We consider a pancake-shaped system in a harmonic 
potential 

Vtrap(r) = ^muilx^ + ^niLoly^ + ^mulz^, 

where the trapping frequencies are LUr = ^ LOy and 
ojz, with ujz ^ uJr- For a sufficiently strong trapping po- 
tential in the ^-direction, the condensate wave function 
and the thermodynamically relevant quasiparticle ampli- 
tudes can be approximated to be in cylindrical coordi- 
nates (r, 9, z) of the factorized form 

Co(r) = Co(r)a(z)e^'"^ (23) 
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and 

u{r) = ^^(r)c^(z)e'(*+")^ (24) 
v{r) = u(r)CT(^)e*(-«^+'")^ (25) 

where a{z) = e^^^/(^°^)/\/ azTT^/^ is a Gaussian profile 
and tti = y^h/mui are the harmonic oscillator lengths 
of the trap. In the following, we consider only the case 
m = of an irrotational condensate. Using Eqs. (23), 
(24) and (25), the Gross-Pitaevskii and Bogoliubov equa- 
tions reduce to equations of the radial coordinate only. 
The z-dcpcndcncc of the GP, the GGP and the Bogoli- 
ubov equations is reduced explicitly by multiplying them 
with a{z) and integrating over z. This results in equa- 
tions similar to the original ones except that the interac- 
tion strength Uq = 4Trh^a/m is replaced with its quasi 
two-dimensional version Uq^ = 2y/2^fuVzaza and the 
chemical potential is shifted by hujz/2. We use a^, l/cOr 
and hujr as units of length, time and energy, respectively. 
In these units the dimensionless interaction strength be- 
comes Uq^ = 2y/2TTa/az- A peculiarity of the reduced 
equations is that they are independent of the trapping 
frequency iVr, and hence our results apply for all cUr, pro- 
vided that Wz ^ U)r- 

In numerical calculations, we use spatial discretization 
and finitc-diffcrcncc methods. The groimd state solu- 
tions of the non-linear GP and GGP equations are found 
by a norm-conserving imaginary time integration method 
based on the Crank-Nicholson scheme. The computa- 
tion speed is enhanced by using a multigrid method, in 
which the grid is made gradually denser during the com- 
putation. On the other hand, in spatial discretization 
the Bogoliubov eigenvalue equation becomes a matrix 
eigenvalue equation, with the coefficic;nt matrix having 
a narrow band. The eigenvalues of this matrix are found 
by implicitly restarted Arnoldi method implemented in 
the numerical library ARPACK. The quasiparticle ampli- 
tudes are calculated from the Bogoliubov equations up to 
some cut-off energy Ecut, above which a semiclassical ap- 
proximation [22] is used for calculating their contribution 
to mean-field potentials. The use of the semiclassical ap- 
proximation enhances the convergence of the results as a 
function of i^cut- 

The anomalous average and the energy shifts AE^ con- 
tain ultraviolet divergences due to the interaction con- 
tact potential approximation, and they have to be renor- 
malized. The renormalization scheme presented in the 
Appendix shows that the divergent part of the anoma- 
lous average is proportional to Co(i")- Using this infor- 
mation we determine the value of the interaction correc- 
tion AUo{Ecut) due to the excitations below Ecut- This 
coefficient also suffices to determine the proper renor- 
malization subtraction for the energy shift AE^^. As an 
approximation to this proper renormalization, the ultra- 
violet divergences can also be removed by neglecting the 
zcro-tcmpcraturc parts of the terms fi;(r) and AE^. For 
dilute gases this should be a good approximation [10]. 
This simpler renormalization method is also computa- 



tionally much faster, since the summations in the en- 
ergy correction AE^ converge and many terms may be 
neglected within computational accuracy. In the more 
accurate renormalization scheme that we have used, the 
summations diverge and hence all the terms up to Scut 
have to be taken into account. 

In calculating the spectral distributions Fp{Lo — ij), we 
have to use a finite imaginary part 7 to avoid divergences 
in Eq. (A. 6). The value of 7 is estimated for each excita- 
tion and temperature separately to be small enough for 
not to affect the mean value of the spectral distribution 
nor its Fourier transform in the regime we have presented 
it. In a finite system, the spectral distributions consist 
of discrete Lorentzian peaks with widths proportional 
to 7. Thus, when calculating the spectral density, the 
smaller the regulator 7 is, the finer grid for the real part 
u must be used. Together with the double summation in 
Eq. (A. 6), this can increase the computational cost of the 
spectrum to be orders of magnitude larger than the cost 
needed in solving the excitations from the Bogoliubov 
equations. To make the computation of the spectrum 
more efficient, the terms Aijk and Bijk in Eq. (A. 6) are 
calculated only once and stored in memory. In addition, 
a comparable speedup is achieved by regrouping the sum- 
mation terms according to their behaviour as functions 
of uj: for slowly varying terms, one can use much sparser 
w-discretization. 

We note that in the end we take the limit 7 0+, 
as was suggested in Ref. [10]. However, in Refs. [2, 16, 
17] a finite value of the order 7 ~ 10~^a;r- was used, 
motivated by the finite experimental observation time, 
and in Ref. [18] the value of 7 was taken to be 5 x lO^'^^Wr- 
However, the inclusion of a finite 7 is only one way to 
model the finite observation time. Since the finite value 
of the regulator 7 has also an unphysical effect of shifting 
the excitation energies, it might not be the best way to 
model the restricted observation time. 



IV. EXCITATION SPECTRA 

In this section, we present and analyze the results for 
the mean energies of excitations as functions of temper- 
ature. We model a pancake-shaped cloud consisting of 
N = 2000 ^■'Na atoms trapped in the tight direction with 
the trapping frequency uiz = 2itx 350 Hz. As pointed out 
in Sec. Ill, the radial trapping frequency 
may be chosen freely with only the constraint cj^ ^ cj,.. 
These parameters are chosen for convenience to coincide 
with the ones used in Ref. [23], in which the excitation 
energies were calculated from the self-consistent HFB- 
Popov theory. Previously, energies of a few modes for 
spherically symmetric systems [17] and condensates con- 
taining a vortex line [18] have been computed within this 
second order theory. 

Figure 1 shows the zcroth order energies of the lowest 
energy excitations at zero temperature. Within the accu- 
racy of the figure, the zeroth order Bogoliubov energies 
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FIG. 1: Excitation energies of the lowest part of the spectrum 
at zero temperature. The energies of excitations marked with 
squares and diamonds are presented as functions of tempera- 
ture in Figs. 2 and 3, respectively. The inset shows the con- 
densate fraction as function of temperature (solid line), and 
the exact result 1 — (T/Tc)^ for the non-interacting system 
(dashdot line). 
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FIG. 2: Temperature dependence of the mean energies of the 
modes marked with squares in Fig. 1. Dots correspond to the 
second order theory, and solid lines to the HFB-Popov theory. 
The dashed line indicates the exact energy huir of the Kohn 
modes. 



coincide with those of the full second order theory. Since 
the condensate is irrotational, the spectrum is symmetric 
with respect to inversion qg —^ —qe, and thus only exci- 
tations with non-negative angular momenta are shown. 
The inset presents the condensate fraction as a function 
of the temperature (solid line) compared with the non- 
interacting gas result N^{T)/N = 1 - {T/T^f (dashdot 
line). We identify the condensation temperature Tc as 
the point where the condensate fraction obtains its max- 
imum second derivative with respect to the temperature. 
The theory is probably not reliable above or in the vicin- 
ity Tc, although we present its predictions also in this 



regime. 



The energies of the modes marked with squares and 
diamonds in Fig. 1 are presented as functions of tempera- 
ture in Figs. 2 and 3, respectively. The mean values of the 
spectral distributions of the excitations within the second 
order theory are shown, in addition to the corresponding 
HFB-Popov results that were obtained by neglecting the 
anomalous average, the second order corrections and the 
terms in the second line of Eq. (A. 3). For the low- lying 
modes, some of the second order energies are observed to 
cross the Popov results as seen in Fig. 2. On the other 
hand. Fig. 3 shows that the second order theory for the 
higher lying modes yields systematically lower excitation 
energies than the Popov theory. The Popov results seen 
in Fig. 2 are consistent with the energies calculated in 
Ref. [23]. 



In Fig. 2, the dashed line corresponds to the energy 
TiLDr of the exact center of mass oscillation modes, the 
Kohn modes. According to the generalization [24, 25] of 
Kohn's theorem [26], a system of harmonically trapped 
interacting particles in any eigenstate of the Hamiltonian 
has an eigenstate with the amount Twji higher energy, i.e., 
the exact diagonalization of the Hamiltonian should yield 
a spectrum with the eigenenergy fiuji. The Bogoliubov 
theory, in which the thermal gas component is neglected, 
implies Kohn modes to have this exact energy. In the 
higher order theories, the dynamics of the thermal gas 
and its interaction with the condensate have to be taken 
into account accurately to obtain results in agreement 
with the Kohn theorem. Figure 2 shows that within the 
second order theory the energy of the Kohn mode is very 
close to fiLOr for temperatures T < 0.8Tc. Taking into ac- 
count perturbation theory terms beyond the second order 
should yield energies even closer to the exact result. 



The lowest mode with vanishing angular momentum 
is the breathing mode corresponding to uniform scal- 
ing oscillations of the condensate. In the case of a two- 
dimensional harmonically trapped gas interacting via the 
contact potential, it has been shown using the scaling 
symmetry of the Hamiltonian that there exists a state 
that has energy 2hujr in excess to the ground state [27]. 
This excitation is identified with the breathing mode. 
The Bogoliubov theory yields exactly the energy 2fiuJr 
for the breathing mode, while the Popov and the second 
order theories do not, as can be seen in Fig. 2. Since the 
interaction potential has to be renormalized and hence 
deviates from that used in Ref. [27] for modes with high 
energy, the applicability of the exact result is somewhat 
questionable at high temperatures, where the physics is 
not determined by the low-lying modes alone. It is shown 
in Fig. 2 that the energy of the breathing mode is lower 
than 2hujr and the deviation from 2hujr increases with 
the temperature. 
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FIG. 3: Temperature dependence of the mean energies of the 
modes marked with diamonds in Fig. 1 within second order 
theory (dots) compared with the HFB-Popov results (sohd 
hues). The energies in the lower block correspond to excita- 
tions with vanishing angular momentum and the upper block 
to excitations with qe « 40. 
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V. DECAY OF THE EXCITATIONS 

In experiments and theoretical studies, the decay of an 
excitation is commonly characterized only by the damp- 
ing rate related to the exponential decay of the oscilla- 
tion amplitude. For infinite systems, the excitation spec- 
trum is continuous and the spectral distributions Fp{ijj) 
of the excitations are Lorentzian peaks, implying indeed 
an exponential decay of the mode oscillations. The mean 
value of the Lorentzian gives the mode frequency and 
its width the damping rate. However, for trapped, finite 
systems the spectrum is discrete and the spectral distri- 
butions generally have more complicated forms. Espe- 
cially, the dynamics implied by these distributions can 
be more complicated than just the simple exponential 
decay. From the computed spectral distributions of the 
oscillations, we have studied the validity of the expo- 
nential decay approximation for the finite system under 
question. 

In fact, the computed spectral distributions consist of 
discrete peaks, as noted also previously [17, 18], and the 
form of the distributions is typically far from a simple 
Lorentzian profile. This seems to imply the excitation 
amplitudes to have a complicated modulation in time. 
The extreme case of this modulation is the collapse and 
revival of the corresponding excitation amplitude. This 
phenomenon can be seen for the breathing mode at zero 
temperature. Figure 4(a) displays the computed spectral 
distribution for the breathing mode, which consists of two 
large, well separated peaks. This implies strong beating 
behaviour in the mode amplitude, seen in Fig. 4(b), in 
which the amplitude of the oscillation collapses in time 
t = Ab/ujr, but revives as time elapses. In fact, the am- 



FIG. 4: (a) Spectral distribution Fiuj) of the breathing mode 
(the lowest mode with qg — 0) at zero temperature and (b) 
its Fourier transform. 



plitude has a beating behaviour with a base frequency 
given by the mean value of the peaks in Fig. 4(a), and a 
beating frequency inversely proportional to the distance 
between the peaks. 

The two peaks in Fig. 4(a) are due to a Beliaev pro- 
cess [28] resonance, in which a breathing mode quasipar- 
ticle with Bogoliubov energy 2huJr decays into two Kohn 
mode quanta with opposite angular momenta and energy 
tkUr. Owing to the temperature independent terms in 
Eq. (A. 6), the Beliaev process may take place even at zero 
temperature. For the simpler, approximate renormaliza- 
tion scheme in which the temperature independent terms 
in the anomalous average and the energy correction AE^ 
are completely neglected, this Beliaev process cannot oc- 
cur at zero temperature, and the spectral distribution 
consists of only a single peak. For the Kohn mode, there 
are no modes into which it could decay via Beliaev pro- 
cesses and hence the oscillation amplitude of the center 
of mass is constant in time at zero temperature. At finite 
temperatures, the resonant Landau process, the inverse 
of the Beliaev process, splits the spectral distribution of 
the Kohn modes. 

The non-trivial spectral distributions of the breathing 
mode and the Kohn modes seem to contradict the results 
for the exact energies of these modes discussed above. 
This is partly due to the accidental strong resonant Beli- 
aev and Landau processes which probably weaken the ac- 
curacy of perturbation theory. If the processes, in which 
the quasiparticles decay into Kohn or breathing modes. 
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FIG. 5: (a) Spectral distribution F{u) of the second lowest 
mode with ge = 1 and (b) its Fourier transform at the tem- 
perature T = 0.64rc. 



are neglected by hand as in Ref. [29], the effects of the 
resonances are removed and the spectral distributions of 
these modes become narrower. However, it is not evident 
that this procedure yields more accurate mean energies. 
In conclusion, a pure collapse and revival of the breath- 
ing mode is probably only an artefact of the second order 
theory — it would be interesting to investigate whether 
higher order calculations would improve the situation in 
this respect. 

The spectral distribution and the dynamics of the 
second lowest mode with qg — 1 a,t the temperature 
T = 0.64Tc is shown in Fig. 5. The distribution is ob- 
viously far from Lorentzian form, consisting of several 
asymmetrically separated peaks, and the dynamics is 
more complicated than the zero temperature result for 
the breathing mode presented in Fig. 4(b). The col- 
lapse and revival behaviour is clearly seen, although it 
is weaker than in Fig. 4(b). Our calculations have also 
showed collapse and ravival of many other elementary 
excitations. 

In Ref. [30], the Beliaev decay has been reported to 
be observed in the case of the scissors mode which cor- 
responds to scissors-like density fluctuation of the con- 
densate. In the experiment, the trapping frequency ratio 
lOz/lOx was adjusted such that the energy i?xz of the scis- 
sors mode in the xz-plane was twice the energy E^y in 
the xy-plane. The amplitude of the oscillation is shown 
in Fig. 3(c) of Ref. [30], and it was observed that the am- 



plitude of the mode decreases and increases in time. It 
was also observed that the strength of this phenomenon 
was peaked into the position of the Beliaev resonance as 
a function of the trap asymmetry ratio ujx/ujy, implying 
that the Beliaev process between these modes is respon- 
sible for this effect. 

Within the second order formalism, one can interpret 
this phenomenon in the following way: Provided that 
the numerator in the term corresponding to the Beliaev 
decay of the scissors mode does not vanish in Eq. (A. 6) 
and other processes are not important, the second order 
theory yields a spectral distribution proportional to 

F,e(w) - lm[l/{nu-Qi-E^., + A/{hLu~Qi-E^, + E,f,))], 

where A is the amplitude of the second order correction 
and -Eofr is the energy indicating how much the Beliaev 
process is off-resonant. If £^ofr = the distribution i^sc(^^) 
consists of two peaks whose distance is determined by 
and the dynamics of the mode corresponds a pure col- 
lapse and revival as shown in Fig. 4 in the case of the 
breathing mode. With increasing i?ofr, one of the two 
peaks becomes smaller and the peaks are shifted in a 
such way that their mean value is kept at i?xz • This cor- 
responds to an oscillation in which the amplitude does 
not vanish completely at any moment, and ultimately 
when i?ofr ^ oo it remains constant, which qualitatively 
explains the observation in the experiments. It is possible 
that due the parity of the scissors modes the second order 
amplitude vanishes. However, the higher order correc- 
tions may still have non- vanishing amplitudes, resulting 
in qualitatively same kind of effect. 



VI. CONCLUSIONS 

We used the gapless second order theory developed in 
Refs. [9, 10] to calculate the excitation energies and dy- 
namics of the collective excitations for a partially con- 
densed, harmonically trapped quasi 2D bosonic gas. The 
results satisfy the Kohn theorem quite accurately for 
temperatures T < O.STc- The energies of the HFB-Popov 
and the second order theory crossed as functions of tem- 
perature for some of the low-lying modes, while the sec- 
ond order theory systematically yields smaller energies 
for the higher lying modes. The first experimental ob- 
servations of the Beliaev damping were discussed within 
the second order theory and it was found that this theory 
qualitatively accounts for the observations. 

The computed spectral densities also imply collapse 
and revival of many elementary excitations. The zero 
temperature spectral distribution of the breathing mode 
is characterized by two large, well-separated peaks and 
the oscillation amplitude consequently displays strong 
collapse and revival behaviour within the second order 
theory. This is due to the resonant Beliaev process, in 
which one quasiparticle in the breathing mode decays 
into two quasiparticles in the Kohn modes with opposite 
angular momenta. The result seems to contradict exact 
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analytical results for the breathing mode energy [27], and 
is probably due to weak convergence of the perturbation 
theory for this mode. The calculations can in princi- 
ple be extended to higher order, but this soon results in 
overwhelming computational difficulties. It would also 
be interesting to upgrade the calculations to be self- 
consistent [9, 10], such that the perturbative energy cor- 
rections are inserted into the eigenvalue equations, which 
are then solved iteratively. 
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APPENDIX: ULTRAVIOLET 
RENORMALIZATION 



Since the low-energy effective contact potential ap- 
proximation for the interactions between the particles is 
not valid at high energies, the bare anomalous average 
K(r) is ultraviolet divergent. To remove this divergence in 
a proper way, following Ref. [16], the interaction strength 
Uo in GGP equation (14) must be replaced by Uq + AUo, 
where 



AUo = Ul 



dk 



(27r)3 n'k^' 



(A.l) 



The term AC/q is divergent, and cancels the divergence of 
the anomalous average. Combining the divergent terms, 
one obtains the renormalized anomalous average 

„^{r)=n{v)+N,^C,l{v), (A.2) 

which is finite. The GGP equation is now properly renor- 
malized provided that the anomalous average is replaced 
by the renormalized one. 

The perturbation Hamiltonian given in Eq. (17) 
is used to calculate the total energy of the system 
E{No, ni, 712, . . . ) up to second order in perturbation the- 
ory for the given quasiparticle distribution {n,}. The ex- 
citation energies are obtained as energy differences Ep = 
E{No - ANp, ni,n2,...,np + l,...)- E{No, m, n2, . . . ), 
where ANp = J dr[|ui(r)p -|- |i'i(r)|^] is the amount of 
particles transferred to the mode p. The energy correc- 



tion terms appearing in Eq. (20) are given by 
AEl = Uo J dr[2p(r)Apj,(r) 

^{2Ap^(r 



-hRe[K*(r)AK(r)] - 
-FRe[AK*(r)AK(r)] } 
Ai^rhape = Wo/dr 

xRe{[Co'(r)-Co'(r)] AK,(r)} 



(A.3) 



+ 2NoUo I dr 
+ 

X / dr 



Co^(r)-Co^(r)] App{r) 
2NoUo J drCo'(r) [up{r) + Vp{r)] 

Co{r)-Ur)\[up{r)+Vp{r)], (A.4) 



AEl = (A - Ag) J drApp(r), (A.5) 



A^f(.') = -E' 



2 1 Apij I ^ 



z' + ei + Ci 



+ 



'^\Bpij\'^ "' 



{l-\-ni+ Uj) 



[ni-nj) 



£7 



-I- 



^'[ffl-^(2n, + l) 



where in the primed summations one excludes the terms 

in which all the summation indices are equal to p. These 
diagonal terms are negligible in the current calculations. 
The contribution to the density of the thermal gas and 
the anomalous average due to the mode p are defined 
as App{v) = \vp{v)\^ + |wp(r)p and AKp(r) = itp(r)w*(r). 
Moreover, we have replaced the energy with a complex 
variable z' . The second order matrix elements are written 
as 

A^jk = Vn^Uo y"dr|Co(r)[wi(r)Mi(r)ufe(r) 

+ Ui {Y)vj {r)uk (r) + Ui {r)uj {r)vk (r)] 
+ Cd(i")K(i")«j(i")"fc(r) 

+ Vi{r)uj{r)vkir) + Vi{r)vj{r)uk(r)]y (A.7) 

Bijk = Vn'oUo j dr{Co*(r)[<(r)u, (r)ufe(r) 

+ <(r)^^j(r)wfe(r) + v*{T)uj{v)vk{T)\ 
+ Co(r)[<(r)Uj(r)t;fe(r) 

+ <W^^j(r)wfc(r) + v*{v)vi{v)vk{v)\ }. (A.8) 
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The energy correction AE^ in Eq. (A. 6) is ultraviolet 
divergent. However, the renormalization of the anoma- 
lous average implies that the bare second order correction 
E^{z) is to be replaced by the renormalized one 

AEP^{z') = El{z') + 2NqAUo 

X |dr|Co(r)|^(|«p(r)|2+K(r)|2).(A.9) 
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